Signatures of Wigner molecule formation in interacting Dirac fermion quantum dots 
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We study N interacting massless Dirac fermions confined in a two-dimensional quantum dot. 
Physical realizations of this problem include a graphene monolayer and the surface state of a strong 
topological insulator. We consider both a magnetic confinement and an infinite mass confinement. 
The ground state energy is computed as a function of the effective interaction parameter a from 
the Hartree-Fock approximation and, alternatively, by employing the Miiller exchange functional. 
For N = 2, we compare those approximations to exact diagonalization results. The Hartree-Fock 
energies are highly accurate for the most relevant interaction range a < 2, but the Miiller functional 
leads to an unphysical instability when a > 0.756. Up to 20 particles were studied using Hartree- 
Fock calculations. Wigner molecule formation was observed for strong but realistic interactions, 
accompanied by a rich peak structure in the addition energy spectrum. 

PACS numbers: 03.65.Pm, 73.22.Pr, 71.15.Bf, 73.21.La 
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I. INTRODUCTION 

Massless two-dimensional (2D) Dirac fermions are of 
central importance in several condensed matter appli- 
cations of current interest, in particular for monolayer 
grapheneiiSi and for the surface state of a 3D strong 
topological insulator (TT)jL&. These systems offer read- 
ily accessible table-top realizations of relativistic quan- 
tum physics, where electron-electron interactions are typ- 
ically much stronger than in atomic physics. Inter- 
actions are characterized by an effective fine structure 
constant a, where for graphene a ~ 1 to 2 depend- 
ing on experimental details^ while in TIs a is proba- 
bly somewhat smaller due to the large dielectric con- 
stant of the relevant thermoelectric materials (for in- 
stance, Bi 2 Se3 or Bi 2 Te3)£ We here study the problem 
of TV massless 2D Dirac quasi-particlcs confined to a cir- 
cular quantum dot of radius R and interacting through 
the Coulomb potential (with prefactor oc a). Quantum 
dots formed in 2D semiconductor heterostructures have 
been studied in much detail over the past two decades, 
both experimentally^ and theoretically^ Given the ex- 
ceptional properties of Dirac fermions and the unique 
properties of the underlying materials, it is of consid- 
erable practical and fundamental interest to investigate 
Dirac fermion quantum dots. Since the commonly em- 
ployed electrostatic gating^ is problematic due to the (re- 
cently observed^) Klein tunneling phenomenon, the ques- 
tion of how to confine Dirac fermions in a controlled man- 
ner arises. While quasi-bound states induced by electro- 
static potentials have also been studied^— we here con- 
sider two types of stable confinement: (i) an infinite- mass 
boundary conditioner— on the single-particle wavefunc- 
tion at r = R, and (ii) confinement by a spatially in- 
homogeneous magnetic field profile*^— Graphene dots 
have already been investigated experimentally by several 
groups,— ~— where confinement was so far created litho- 
graphically. While this (approximately) corresponds to 
case (i) above, such a procedure may give rise to uncon- 



trolled disorder effects along the boundary, and route (ii) 
may offer a promising alternative for future experiments, 
see also Ref. H3- For the TI surface state, we arc not 
aware of experimental reports of quantum dot physics, 
but confinement should be achievable as well using, e.g., 
suitably arranged close-by ferromagnetic layers. 

On the theoretical side, another difficulty arises from 
the Dirac nature of the quasi-particles when one at- 
tempts to include electron-electron interactions. For 
the TV-particle problem, where a first-quantized formu- 
lation generally offers the most natural routed the prob- 
lem arises from the unboundedness of the single-particle 
Dirac Hamiltonian in Eq. ((T|) below. This causes the 
Brown-Ravenhall 2 ^ "disease:" the unbounded spectrum 
allows particles to lose arbitrary amounts of energy by 
transferring their energy in (real) scattering events to 
other particles. To circumvent this problem, suppose 
that the chemical potential is located just above the Dirac 
point. We then follow Sucher— and confine the Hilbert 
space to positive-energy eigenstates of the full single- 
particle problem, i.e., we assume an inert filled Dirac sea. 
This projection approach has been successfully employed 
in the same context before ; 19 ' 20 and one can analyze also 
other values for the chemical potential. The accuracy 
of this method was carefully assessed in Ref. [l^. In 
short, the presence of a spectral gap due to confinement 
allows to implement Sucher's no-pair approximation 2 ^ 
since electron-hole pair excitations neglected in this ap- 
proach have to overcome the gap. 

Below we show and compare results from three differ- 
ent computational approaches. In particular, we perform 
self-consistent Hartree-Fock (HF) calculations and, in ad- 
dition, study a similar self-consistent variational proce- 
dure using the so-called Miiller density matrix functional 
(replacing the Fock term) i 30 i 31 In atomic physics applica- 
tions, the Miiller functional is sometimes superior to the 
HF approach and is valuable because it yields a lower 
bound for the ground-state energy^ We note that HF 
calculations for graphene dots with infinite-mass confine- 
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merit have also been carried out by other groups! 32 ' 33 
While the HF approximation is known to provide an 
upper bound for the exact (within Sucher's projection 
approach) ground-state energy^ the Miiller functional 
is again expected to generate a variational lower bound. 
We will compare results from those two approaches to ex- 
act diagonalization computations for N = 2 interacting 
Dirac quasi-particlcs. The model and those three numer- 
ical approaches are described in Sec. UH while the com- 
parison for N = 2 can be found in Sec. IIIII As expected, 
the exact results are always bracketed by the results ob- 
tained from the Miiller functional and under the HF ap- 
proximation. However, within the range a < 2 studied in 
this work, we find that the HF results are much closer to 
the exact results and provide a rather accurate approxi- 
mation. However, results based on the Miiller functional 
show an unphysical divergence when a > 0.756, and are 
less accurate than the HF results for small a. (Of course, 
for a — > 0, all three methods recover the correct nonin- 
teracting results.) Having established the HF approach 
as highly accurate approach for a < 2 and N — 2, we 
continue in Sec. HVI with a presentation of HF results for 
N > 2 Dirac particles in a quantum dot with infinite- 
mass confinement. Besides the ground-state energy, we 
study various physical obscrvables like the particle den- 
sity or the spin density. Our results suggest that in a 
confined geometry Dirac particles can form a "Wigner 
molecule" as previously discussed for Schrodinger par- 
ticles in semiconductor dots.— i 3 ^— When the Coulomb 
interactions dominate over the kinetic energy, a Wigner 
crystal can be formed where electrons spontaneously or- 
der in a crystalline structure. The presence of a con- 
fining potential makes this Wigner crystallization more 
favorable,— and although no Wigner crystal is expected 
for bulk graphenej22 we find that the confined geome- 
try allows for a finite-size Wigner "molecule" even for 
Dirac fermions. The paper concludes with a discussion 
in Sec. EJ 

II. MODEL AND COMPUTATIONAL 
APPROACHES 

In this section we discuss the model studied in this 
work, and address the different calculational schemes em- 
ployed to study the iV-particle problem for interacting 
Dirac fermions in a quantum dot. 

A. Single-particle model 

We consider a single species of massless 2D Dirac 
fermions described by the single-particle Hamiltonian 
(— e < is the electron charge) 

H = v F cr- (p+-A) +Ma 3 , (1) 

where er = (01,02) an d the Pauli matrices Oi refer to 
the sublattice structure of the honeycomb lattice for 



graphcno^i or to the electronic spin degree of freedom for 
the TI surface stated The Fermi velocity in graphene is 
vf ~ 10 6 m/s, while the corresponding value for the 
TI surface state is approximately half this value. A 
single Dirac cone as in Eq. (Q]) can be realized for a 
TI surface^ but in graphene there generally is a four- 
fold degeneracy due to the valley and spin degrees of 
freedom^ For graphene, we then assume a spin- and 
valley-polarized situation where the single-valley theory 
[Eq. (J!}] gives useful predictions^ In fact, our basic qual- 
itative conclusion, i.e., Wigner crystallization is possible 
in graphene dots, is also found from HF calculations in- 
cluding the spin and valley degrees of freedom^ In ad- 
dition, we allow for a static vector potential A(r) cor- 
responding to inhomogeneous magnetic fields or, in the 
case of graphene, also to strain-induced pseudo-magnetic 
fields.— Finally, M(r) corresponds to a mass term. In 
order to form a quantum dot, where Dirac fermions are 
confined to a bounded spatial region, say, a disk of radius 
R around the origin, we now consider the two possibil- 
ities mentioned in Sec. HI We study circularly symmet- 
ric cases, where the total angular momentum operator 
J = —ihdcf, + ha z /2 is conserved and has eigenvalues Hj 
with half-integer j = m + 1/2, m G Z. Single- particle 
solutions to HqiP = Etjj can then be written as 

^c.«=^(A"Lm)' (2) 

In what follows, we measure energies (lengths) in units 
of Hvf/R {R), and we always focus on E > solutions. 

(i) Infinite mass confinement. — A well-known way 
to describe confinement theoretically is to impose an 
infinite-mass boundary condition on the wavefunction 
[Eq. ©] atr = 1, i.e., M(r < 1) = and M(r > 1) -> 00. 
As shown by Berry and Mondragon^i the effect of M(r) 
in Eq. ([I} is then fully captured by the boundary condi- 
tion 

^l,m(l) = V>2,m(l), (3) 

stating that no current flows through the boundary. With 
the Bessel function J m , the Dirac equation is solved for 
r < 1 by the Ansatz 

ipi,m(r) = AJ m (Er), ip2,m(r) = AJ m+1 (Er), 

where the boundary condition ([3]) yields the energy quan- 
tization condition^ 

■^m(^mn) ^m+l(^m)i)- (4) 

This equation has to be solved numerically. (Note again 
that E is given in units of Hvf/R and r in units of R.) 
Positive-energy solutions for given m are then labeled by 
?i € N. We mention in passing that there are no zero- 
energy solutions^ The normalization factor A is 

A m n = [ 7r (^m — Jm— 1 Jm+1 + J'm+l ~ JrnJm+2 )] ^ 2 j (5) 
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where all Bcsscl functions are evaluated at E mn . To sum- 
marize, the single-particle solutions ip a under a circular 
infinite-mass confinement are labeled by a — (m, n) with 
in E Z and n E N. The eigenenergies E a follow by solving 
Eq. ((4]) and the eigenspinor is (r < 1) 

For a detailed discussion of the single-particle spectrum, 
see Ref. GJ. 

(ii) Magnetic confinement. — A second possibility to 
confine Dirac quasi-particlcs is to employ spatially in- 
homogeneous magnetic fields. This possibility has been 
explored theoretically before,— ~— and we study the sim- 
plest case of a piecewise constant magnetic field, B(r) = 
BQ(r — 1) with B > and the Heaviside step function 
Q(x). The eigenenergies for this single-particle problem 
can be found numerically and were given in Ref. |2CI The 
spectrum contains "dot states" , with probability den- 
sity concentrated in the disk region r < 1, plus rcl- 
ativistic bulk Landau states for r > 1. The Landau 
states are weakly perturbed by the presence of the di- 
mensionless "missing flux" parameter S := R 2 /2£ 2 , where 
t := \JcjeB is the magnetic length. Because of this 
perturbation, Landau level energies are slightly shifted 
away from their standard bulk value, but dot states can 
be clearly distinguished in the single-particle spectrum. 
With chemical potential chosen such that all bulk Lan- 
dau states below the first Landau state, E^> := \/2R/i 



The coefficient Ck^i vanishes when I + |fc| is odd or when 
I < \k\. For k = I = 0, we have Ck.i = 1/2. In all 
remaining cases, we obtain 

_ (2/-l)!! ( ™ /2 {n -l/2)(n-l-l) 
W 2'+!?! J- -I- n(n-l- 1/2) ' 

n— 1 v ' 1 

Equation (jSJ is then evaluated by numerical integration 
routines and yields the interaction matrix elements. 

B. Numerical approaches 

Next, we briefly describe three different numerical ap- 
proaches to obtain the ground-state energy for a quan- 
tum dot containing TV Dirac fcrmions, namely HF simula- 
tions, the Miiller density matrix formulation, and exact 



(in units of hvp /R) , are filled, the relevant dot states are 
in the window < E a < E^K All cigenstates can again 
be labeled by a — (m,n), i.e., using angular momentum 
j = m+ 1/2 and the index n £ N. For given missing flux 
5, there are N b (6) dot states, where N b increases with 
increasing S, see Ref. [2(1 The TV-particle problem can 
then be studied for N < Nb(S) only. In fact, due to the 
repulsive interactions, the maximum number of bound 
electrons may be lowered even further For the infinite- 
mass confinement [case (i)], there is no constraint on the 
number of particles held by the dot. 

We now add electron-electron interactions to the 7V- 
particle problem. The Coulomb interaction matrix ele- 
ments are given in terms of the cigenspinors ij) a , 

Vaa'b'b :=aj (V4 • A) (r) {f a , ■ Vv) (O, (7) 

with V a 'abb' = Vaa'Vb- Due to total angular momentum 
conservation, only matrix elements with m a + m a i = 
m + my do not vanish. Interaction matrix elements 
with large momentum exchange k = mb ~ m a (k € Z) 
are numerically small^ but all possible values of k (for a 
chosen basis size) are taken into account below. For the 
magnetic dot [case (ii)] , the matrix elements ([7]) are most 
conveniently evaluated by expanding ip a in conventional 
relativistic Landau level states^ For the infinite-mass 
confinement [case (i)] , after inserting Eq. © into Eq. Q , 
some algebra [cf . also Appendix B of Ref. [2(| yields 
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diagonalization (for N = 2). In order to have a well- 
defined many-body problem, we follow Sucher— and re- 
strict ourselves to the projected single-particle space, i.e., 
we assume an inert filled Dirac sea. Hence summations 
over a = (m, n) will only include positive-energy single- 
particle solutions (E a > 0). In the numerical calcula- 
tions, the basis size (i.e., the number K of single-particle 
orbitals spanning the Hilbcrt space) was always chosen 
sufficiently large to ensure convergence. Failure to con- 
verge indicates an instability of the method, as we will 
see in the case of the Miiller functional for strong inter- 
actions. 



First, the Hartrcc-Fock approach amounts to the self- 
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consistent minimization of the functional 
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(Vaa'b'b ~ Vaa> bb')l a> b' 1 ab , 



(9) 

where the density matrix 7 obeys 7 2 = 7 and tr(7) = N . 
In our case, 7 is a real symmetric matrix. The numer- 
ical algorithm to obtain the HF ground state is stan- 
dard and can be found, for instance, in Ref. |2fJ. Second, 
the Miiller density matrix formulation employs a differ- 
ent form for the exchange term, where one minimizes the 
functiona l 30 ' 31 

E M[l] = ^E a "f aa + )^ ^ (Vaa'b'bla'b'Jab (10) 
a aa'bb' 

- Ka' hb <(7 1/2 W(7 1/2 )afc), 

where 7 is again a real symmetric matrix with tr(7) = N , 
but now 7 2 < 7. A stable numerical approach to mini- 
mize £m [7] in Eq. ([TO)), the so-called projected gradient 
algorithm, has been formulated and tested before i 40 ' 41 
We have employed precisely the same method here. Fi- 
nally, the exact numerical diagonalization of the full 
many-body problem is only possible for small particle 
numbers due to the exponential increase in computa- 
tional complexity with increasing N. We have there- 
fore carried out exact diagonalization calculations only 
for N — 2 Dirac fermions, primarily to check the accu- 
racy of the two computationally less expensive but ap- 
proximate alternative approaches. Details of the exact 
diagonalization approach have been described in Ref. fl9L 



III. COMPARISON OF METHODS: N = 2 

In this section, we show and compare the results of the 
three approaches described in Sec. [TT] for N = 2 Dirac 
fermions. For the infinite-mass confinement case, results 
for the ground-state density E(a) are shown in Fig. [JJ 
Clearly, the numerically exact result obtained from ex- 
act diagonalization is bracketed by the HF prediction 
from above and by the Miiller result from below. The 
HF approximation provides very accurate estimates for 
E(a), while the Miiller functional is only reliable for very 
small a. In fact, the application of the Miiller functional 
to the free-space case reveals an intrinsic divergence for 
strong interactions a > a c , where the critical value is 
(see Ref. H^, correcting an earlier attempt 4 ^) 
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0.756, 
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Although this critical value was derived for the case of 
vanishing magnetic field, we anticipate that it applies also 
to the confined geometry, with or without magnetic field, 
since it arises from the fact that both the kinetic energy 
and the Coulomb singularity scale as inverse length for 
short distances, i.e., a regular magnetic field is clearly 
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FIG. 1: (color online) Interaction contribution to the ground- 
state energy, E(a) — E(0) (in units of Twf/R), vs fine structure 
constant a for the infinite-mass confined dot containing N = 
2 particles. The main panel shows the results of the three 
different approaches, see main text. The HF results are very 
close to the exact diagonalization results, while the Miiller 
functional gives a lower bound. Straight lines are a guide to 
the eye only. Inset: Same for Miiller functional, with different 
basis size K. Note the absence of convergence for large a. 
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FIG. 2: (color online) Same as Fig.fTJbut for the magnetically 
confined case with I = R. The basis size K here corresponds 
to 2|m m in|. 



irrelevant for the singular behavior. (Such a result has 
recently been established in a related situation^ 4 -) For 
a > a c , the exchange part in the Miiller functional 
[Eq. ([TO)) ] provides a strong attraction which effectively 
forces particles to form a droplet. In our case, this singu- 
lar behavior implies that the "ground-state" energy drops 
to —00. In numerical computations, this is reflected by 
the fact that the energy becomes cutoff-dcpcndcnt, going 
to —00 as the basis size K grows. This phenomenon is 
clearly visible in the inset of Fig. [1] but a precise com- 
parison of the predicted critical value for a c [Eq. (|11[) ] 
with numerics is difficult. This singularity is an unphys- 
ical artefact of the Miiller density matrix approach and 
indicates that it is only useful for a <C 1. On the other 
hand, the HF approximation is very close to the exact 
value even for a = 2. 




FIG. 3: (color online) HF results for the addition energy A(JV) 
[Eq. (|12[) ] vs particle number N for a dot formed by infinite- 
mass confinement. Results are shown for a — 2 (red circles) 
and for a — (black squares); straight lines are a guide to 
the eye only. Inset: HF results for the energy E(N) vs N, for 
the same interaction parameters. 



A very similar picture emerges from the corresponding 
study of the magnetically confined dot, see Fig. [2] In 
both cases and for all a < 2, the interaction energy ob- 
tained under the HF approximation is less than 1% above 
the corresponding exact value. In the remainder of the 
paper, we will then study N > 2 particles using the HF 
approach. We have compared the results of the Miillcr 
functional for N > 2 to the corresponding HF results as 
well, and with increasing N they come closer. Hence we 
expect that the relative accuracy of the HF results (at 
the least) does not deteriorate for N > 2. 



IV. HARTREE-FOCK RESULTS FOR N > 2 
PARTICLES 

In the previous section, we have established that HF 
calculations are able to provide very accurate estimates 
for the ground-state energy of Dirac fermions in a circular 
quantum dot. In this section, we describe the results of 
our HF calculations for up to N = 20 particles. For 
clarity, we focus on the infinite-mass confinement case, 
but qualitatively similar results were also found for the 
magnetic confinement. 

Figure [3] shows HF results for the TV-dependent addi- 
tion energy, 

A(N) :=E(N + 1)+E(N-1)-2E(N), (12) 

both for a = 2 and for the noninteracting case (a = 0). 
The HF ground-state energy E(N) obtained from our 
self-consistent numerical calculation is shown in the in- 
set of Fig. [3] A peak in the addition energy for some 
N implies a higher stability of the iV-particle dot. In 



FIG. 4: (color online) HF results for the radial density profile 
p(r) vs r for N = 9 particles with several a. 



analogy to atomic and nuclear physics, this N is often 
referred to as "magic number."— While already the non- 
interacting dot has some structure in the addition energy 
spectrum (due to the single-particle spectrum), e.g., the 
small peaks at N = 7 and N = 11 visible in Fig. [3J 
the interacting case is characterized by more pronounced 
features. For a = 2, we observe clear peaks, see Fig. [3j 
corresponding to the magic numbers N = 4,7, 11, 13, 15 
and 18. Although some of these numbers coincide with 
the noninteracting ones, it is evident that the addition en- 
ergy spectrum is drastically changed by electron-electron 
interactions in such a finite-size system. 

The resulting ground-state density p(r) is rotation- 
ally invariant and can therefore be analyzed in terms 
of the angular-averaged density p(r), which is normal- 
ized as f rdrp{r) = N. Figure 2] shows HF results for 
the density p(r) for N — 9 and several a. In the non- 
interacting case (a = 0), the density profile is rather 
smooth, but with increasing a the particles are pushed 
towards the boundary and form a ring. When compar- 
ing the shoulder-like feature apparent in Fig. [4] (around 
r w 0.5) to the corresponding correlation plot (see be- 
low), we find that no significant particle weight is con- 
tained in the shoulder, i.e., with high probability all par- 
ticles are close to the boundary. For N = 19 particles, 
a richer structure emerges, see Fig. [5j with three differ- 
ent spatial "shells" emerging for strong interactions. In 
particular, by integrating over the shown density curve, 
we find that one particle is located near the origin, three 
particles are contained in a second shell around r ~ 0.45, 
and the remaining 15 particles are close to the boundary. 

To obtain more detailed insight we next study the 
density-density correlation function 

g(r,r') = (p(r)p(r')}, (13) 

where r 1 is kept fixed. Monitoring g(r,r') as a function 
of r, the spatial arrangement of the particles in the dot 
can be revealed. 
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FIG. 5: (color online) Same as Fig. g]but for N = 19. 



2D correlation plots for TV = 9 and TV = 19 (with 
a = 2) are shown in Fig. |6] and Fig. [7j respectively. In 
these plots, we keep r' = (0.95, 0) fixed and show the 
correlations as function of r = (x, y) within the dot. For 
TV = 9, Fig. [5] is consistent with all electrons being ar- 
ranged equidistantly on a ring close to the boundary. The 
correlation plot in Fig. [7] for TV = 19 particles also con- 
firms the conclusions reached from the analysis of the 
density plot in Fig. [5l The outermost spatial shell (near 
the boundary) holds 15 particles, a second ring contains 
3 particles, and one particle is located at the center. The 
combined analysis of density and correlation plots for all 
particle numbers under study, TV < 20, results in the shell 
filling sequence in Table |TJ 
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FIG. 6: (color online) Correlation plot g(r, r') for TV = 9 par- 
ticles and a — 2, corresponding to Fig. [4] The position r' is 
fixed at (0.95, 0), and the color scale indicates the correlation 
degree for different r within the quantum dot. 
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FIG. 7: (color online) Same as Fig. H but for TV = 19. 



These observations provide a signature for the onset of 
Wigner molecule behavior, i.e., we have a finite-size sys- 
tem where Wigner crystallization sets in but quantum 
fluctuations are still importantj 35 ' 36 In order to compare 
to the deep Wigner crystallized limit, we now briefly dis- 
cuss the classical limit (which here is defined by taking 
the limit a — > oo), where the electrostatic energy domi- 
nates completely and the kinetic energy can be neglected. 
The repulsive interaction then tries to maximize the dis- 
tance between particles, leading to the formation of spa- 
tial shells. The shell filling sequence for a harmonically 
confined Wigner molecule (of Schrodinger fermions) is 
well known ) 7 ' 34 i 35 and HF calculations have been able 
to capture the Wigner molecule formation^ For the 
2D circular hard-wall confinement considered here, how- 
ever, a different shell filling sequence follows by mini- 
mization of the classical electrostatic energy E C (N) with 
respect to all particle positions r.; = i n within the disk 

(n<R=i), 



E C {N) 



N 

E 



o 



(14) 



A first possible configuration has all particles arranged 
equidistantly on a unit circle, resulting in the classical 
energy 



^(N-l)/2 i 



TV odd, 



AT I V^l ~- 

Z { l^k=l sin(7Tfc/A0 + 2' iV CVCn - 

(15) 

If instead one particle resides at the origin plus TV — 1 
particles on the outer ring as above, the energy of this 
second configuration is 

(TV) = (TV - 1) + (TV - l)a. (16) 



For TV < 16, numerical minimization of Eq. (fT4"]) shows 
that these two configurations always yield the lowest- 
energy solutions. In particular, E^ < E^ for TV < 12, 
see Table HI For 16 < TV < 20, an additional inner ring 
is formed containing TVf 1 > 1 particles, surrounded by 
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TABLE I: Shell filling sequence for a 2D interacting Dirac 
fermion dot with circular hard-wall confinement. Nf de- 
notes the number of particles in the ith spatial shell obtained 
from the minimization of the classical electrostatic energy 
[Eq. (|14|l ]. Ni is the corresponding HF quantity for a = 2, 
see main text. 
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FIG. 8: (color online) HF results for the spin density s r in 
radial direction vs r for N = 9 and various a. 



the angular variable. We find that only the radial com- 
ponent s r (r) := s ■ e r (with e r = r/r) does not vanish. 
The resulting nontrivial spin texture is shown in Fig. [5] 
for N = 9 and several a. 



the outer ring containing iVf = N — N% particles. For 
all N < 20, the classical lowest-energy solution thus has 
at most two spatial shells, but configurations with three 
shells as observed in the quantum calculation are energet- 
ically quite close. The agreement between the shell fill- 
ing sequence observed for a = 2 and in the classical limit 
is not perfect but indicates that we are already rather 
close to the classical limit for a = 2 and have a Wigncr 
molecule, despite of the theoretically predicted absence 
of Wigner crystallization in bulk graphene^i Even for 
a = 1, the above density plots suggest that incipient 
Wigner molecule behavior can be observed. (Of course, 
this is a smooth crossover and not a phase transition.) 
However, the fact that there are still substantial quan- 
tum fluctuations for a = 2 is also clear from the addition 
spectrum in Fig. |3l In the deep classical limit, there is 
much less pronounced structure in the addition energy 
spectrum. 

Finally, we point out that there is also interesting spin 
texture in such a quantum dot. The Pauli matrices in 
Eq. ([1]) are directly connected to the electronic spin den- 
sity in a topological insulator surface via the relation^^ 

s(r) = {s x ,s y ) T = ^(e z x er). (17) 

For the case of graphene, the Pauli matrices refer to the 
sublatticc degree of freedom, which is not easily accessi- 
ble to experiments. The spin density (|17[) points within 
the 2D plane and is always isotropic, i.e., independent of 



V. DISCUSSION 

In this paper, we have discussed interaction effects in 
circular 2D quantum dots where the particles are mass- 
less Dirac fermions. Physical realizations of the studied 
model are given by graphene and the surface state of 
a topological insulator. For the case of two particles, 
we have compared three different methods to establish 
that Hartree-Fock calculations provide highly accurate 
results for physically relevant interaction strengths. An 
alternative method based on the Miiller density matrix 
functional was also studied, but since Muller's Ansatz for 
the two-particle density respects the right normalization 
condition but sacrifices its positivity, it suffers from an 
unphysical divergence for sufficiently strong interactions. 
An improvement would have to take this drawback into 
account, while not dropping the sum rule for the density 
and the convexity of the functional. 

The case of N < 20 particles has then been studied 
using Hartrcc-Fock simulations. The resulting addition 
spectrum of the quantum dot reveals pronounced magic 
numbers that cannot be explained by a noninteracting 
picture. Moreover, the density profiles and the density- 
density correlation functions show that we are rather 
close to the classical limit already for experimentally rele- 
vant interaction parameters (a 1 to 2). The formation 
of spatial shells is a clear signature of a Wigner molecule, 
and we therefore predict that in such a finite-size system 
the usual argument^! for absence of Wigner crystalliza- 
tion of Dirac fermions can be effectively circumvented. 
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